###############################################
#FIGURE 3
###############################################

library(estimatr)
library(ggplot2)
library(grid)

####load cleaned data####
data1 <- read.csv("./data/data1.csv")
data2 <- read.csv("./data/data2.csv")
data3 <- read.csv("./data/data3.csv")


####results without controls####

lm1 <- lm_robust(natid_pca ~ as.factor(trt), data=data1)
lm2 <- lm_robust(natid_pca ~ as.factor(trt), data=data2)
lm3 <- lm_robust(natid_pca ~ as.factor(trt), data=data3)

com1 <- c("Identity Threat\nvs. Control", lm1$coefficients[2], (lm1$coefficients[2] + (1.65*lm1$std.error[2])), (lm1$coefficients[2] - (1.65*lm1$std.error[2])),
          (lm1$coefficients[2] + (1.96*lm1$std.error[2])), (lm1$coefficients[2] - (1.96*lm1$std.error[2]))) 
com2 <- c("Whataboutism\nvs. Identity Threat", lm2$coefficients[2], (lm2$coefficients[2] + (1.65*lm2$std.error[2])), (lm2$coefficients[2] - (1.65*lm2$std.error[2])),
          (lm2$coefficients[2] + (1.96*lm2$std.error[2])), (lm2$coefficients[2] - (1.96*lm2$std.error[2]))) 
com3 <- c("Whataboutism\nvs. Control", lm3$coefficients[2], (lm3$coefficients[2] + (1.65*lm3$std.error[2])), (lm3$coefficients[2] - (1.65*lm3$std.error[2])),
          (lm3$coefficients[2] + (1.96*lm3$std.error[2])), (lm3$coefficients[2] - (1.96*lm3$std.error[2]))) 
com <- data.frame(rbind(com1, com2, com3))
colnames(com) <- c("treatment", "effect", "cih90", "cil90", "cih95", "cil95")
level_order <- c("Identity Threat\nvs. Control", "Whataboutism\nvs. Identity Threat", "Whataboutism\nvs. Control") 

com$effect <- as.numeric(com$effect)
com$cih90 <- as.numeric(com$cih90)
com$cil90 <- as.numeric(com$cil90)
com$cih95 <- as.numeric(com$cih95)
com$cil95 <- as.numeric(com$cil95)
coma <- com


####results with controls####

lm1 <- lm_robust(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1)
lm2 <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)
lm3 <- lm_robust(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data=data3)


com1 <- c("Identity Threat\nvs. Control", lm1$coefficients[2], (lm1$coefficients[2] + (1.65*lm1$std.error[2])), (lm1$coefficients[2] - (1.65*lm1$std.error[2])),
          (lm1$coefficients[2] + (1.96*lm1$std.error[2])), (lm1$coefficients[2] - (1.96*lm1$std.error[2]))) 
com2 <- c("Whataboutism\nvs. Identity Threat", lm2$coefficients[2], (lm2$coefficients[2] + (1.65*lm2$std.error[2])), (lm2$coefficients[2] - (1.65*lm2$std.error[2])),
          (lm2$coefficients[2] + (1.96*lm2$std.error[2])), (lm2$coefficients[2] - (1.96*lm2$std.error[2]))) 
com3 <- c("Whataboutism\nvs. Control", lm3$coefficients[2], (lm3$coefficients[2] + (1.65*lm3$std.error[2])), (lm3$coefficients[2] - (1.65*lm3$std.error[2])),
          (lm3$coefficients[2] + (1.96*lm3$std.error[2])), (lm3$coefficients[2] - (1.96*lm3$std.error[2]))) 
com <- data.frame(rbind(com1, com2, com3))
colnames(com) <- c("treatment", "effect", "cih90", "cil90", "cih95", "cil95")
level_order <- c("Identity Threat\nvs. Control", "Whataboutism\nvs. Identity Threat", "Whataboutism\nvs. Control") 

com$effect <- as.numeric(com$effect)
com$cih90 <- as.numeric(com$cih90)
com$cil90 <- as.numeric(com$cil90)
com$cih95 <- as.numeric(com$cih95)
com$cil95 <- as.numeric(com$cil95)
comb <- com

####combined figure####

coma$Controls <- "Without"
comb$Controls <- "With"
full <- rbind(coma, comb)

level_order <- c("Identity Threat\nvs. Control", "Whataboutism\nvs. Identity Threat", "Whataboutism\nvs. Control") 

ggplot(full, aes(x = factor(treatment, level = level_order), y = effect, color=Controls)) +
  geom_point(stat="identity", size=4, position = position_dodge(width = 0.15)) + 
  geom_linerange(aes(ymin = cil95, ymax = cih95), 
                 position = position_dodge(width = 0.15,),
                 linewidth = 1) +
  geom_linerange(aes(ymin = cil90, ymax = cih90), 
                 position = position_dodge(width = 0.15),
                 linewidth = 2) +
  ylab("Treatment Effect") +
  xlab("Treatment Type") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  theme_bw()+
  theme(
    axis.title=element_text(size=18), 
    axis.title.x = element_text(margin = unit(c(7, 0, 0, 0), "mm")),
    axis.title.y = element_text(margin = unit(c(0, 7, 0, 0), "mm")),
    plot.title = element_text(size=26),
    axis.text=element_text(size=16)
  )  

ggsave("fig3.png", width = 10, height = 6, dpi = 300)

####additional statistics related to these results reported in-text####

#identity threat vs control effect size in SDs
data_clean <- read.csv("./data/data_clean.csv")
datac <- subset(data_clean, data$trt == "control")
lm1$coefficients[2]/sd(datac$natid_pca, na.rm=T)
